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THERMAL CONVECTION AT INFINITE PRANDTL NUMBER 


I. INTRODUCTION 

The problem of determining theoretically the heat transport through a con- 
vecting layer of fluid bounded by horizontal conducting surfaces has received 
much attention lately. Theoretical approaches to this problem may be placed in 
two categories; (1) statistical theories, and (2) direct integration of the fluid's 
equations of motion. Statistical theories seek to determine some of the statisti- 
cally averaged properties of the flow field without computing completely the 
turbulent temperature and velocity fields. Its virtue is that average properties 
are smooth functions of space-time coordinates so that fewer parameters are 
needed for their description than for the chaotic velocity and temperature fields. 
On the other hand, direct integration of the equations of motion is certainly to be 
preferred in those situations in which the flow field is non turbulent, and spacially 
regular; it is also of interest to attempt a solution of the turbulent flow situation 
by this second avenue. A reason for this statement is that current statistical 
theories are not able to estimate accurately their errors. A suitably designed 
integration procedure may thus be used not only to study the flow fields directly, 
but also to assess the accuracy of proposed statistical theories. 

Applications of some of the statistical theories to the thermal convection 
problem has been made for the case in which the Prandtl, number, <r, is put 
equal to infinity (Herring, 1969). The statistical procedures studied included 
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the direct interaction approximation (Kraichnan (1964) ) and the quarinormal 
approximation (an extension of that proposed by Tatsumi, Proudman, and Reid 
(1954), and (1957)). Of these two, only the direct interaction appears to be a 
physically acceptable procedure. It was found to reproduce accurately some 
statistical results constructed from the direct integration approach. The study 
was restricted to moderately small values of Rayleigh number (R £ 10 4 ) and to 
the free boundary conditions. 

An important question raised by the above study is which of the two methods - 
the statistical approximation or a direct integration of the amplitude equations of 
motion is most economical from the computational point of view. For the cr - o o 
limit, it appears that the direct integration method may be at least as feasible as 
the statistical approach. The reason for this statement is that in this limit the 
velocity field does not boundary layer, and consequently the statistical wave 
number cascade of entropy and energy to large horizontal wave numbers is not 
very severe. The velocity and temperature fields may then be represented by 
relatively few horizontal wave number vectors. 

A meaningful amplitude calculation must not only contain a sufficient number 
of mesh points (or Fourier modes) to accurately represent the small scale 
structure of the turbulent flow, but must as well represent its homogeniety and 
isotropy. It appears that present day computors can deal with such a problem 
provided the degree of turbulent activity is not too intense. 
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The present paper presents some results of a direct integration of the 
Boussinesc - Navi er- Stokes equations for the thermal convection problem. The 
Fourier .ransform of the equations of motion are integrated forward in time 
until cjrtain time average properties in the system become steady in time. The 
initial data for the temperature and velocity field are selected from a Gaussian 
d.stribution having horizontal homogeniety and isotropy. Results are obtained 
for the case in which the Prandtl number, cr, is put equal to infinity, which sim- 
plifies the dynamics of the flow by effectively eliminating the vVv nonlinearity 
from the momentum equation. The resulting equations should be valid provided 
the Reynolds number of the flow is small. This condition does not preclude turbu- 
lent flows, since the Pechlet number may be large. Our calculations here is a 
preliminary step toward considering the more complicated case of finite cr con- 
vection. 

The use of Fourier mode techniques permits an entirely energetically con- 
sistent treatment of the nonlinear terms, and, more importantly affords a con- 
venient basis for representing flow fields having dimensions of homogeniety. 

The present problem is set up for a horizontally infinite system, and the initial 
data are consistent with horizontal homogeniety. 

Previous numerical calculations for thermal convection have been made by 
Fromm (1965) and Deardorff (1964). Our calculations differ from theirs in several 
respects. First, the present results for a = co, whereas Fromm and Deardorff put 
o- = 1. The system we study is therefore, in this respect, simpler. However, under 
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these simplified conditions, we carry out a three dimensional integration of the 
equations of motion, where the above authors restrict their results to two dimen- 
sions. Moreover, the method of integration is designed to treat (approximately) 
the flow situation in which the velocity and temperature fields have horizontal 
isotropy and homogeniety. 

The results may be compared to experiments for high Prandtl number 
fluids. Recent experiments for such fluids have been reported by Rossby (1966), 
Somerscales and Dropkin (1966), Dropkin and Somerscales (1965), and most re- 
cently by Somerscales and Gazda (1968). The results of Somerscales and Gazda 
are of particular interest to us in that they present detailed measurements of 
the mean temperature field and (r.m.s.) temperature fluctuation fields for a 
range of Rayleigh numbers (7.39 x 10 s < R < 3.21 x 10 8 ). A curious feature 
of the latter authors results is that the measured mean temperature field is 
asymetric about the midplane of the flow. Such an asymetry is compatible with 
the Boussinesq equations for a horizontally infinite system, but it is difficult to 
rationalize on the basis of theories of turbulent convection. For example, it is 
well known that hexagonal cellular convection produces an asymetric mean tem- 
perature field, but if the flow is time-fluctuating (as indeed the authors report) 
it is difficult to see how the basic hexagonal plan form can persist for a long 
time. Even if the flow is static, the hexagonal plan form may not be stable. 
Asymetrics of the sort reported by Somerscales and Gazda have also been re- 
ported by Somerscales and Dropkin, and by Dropkin, but not by Rossby. 
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The results of the present calculation are in good agreement with the above 
cited experiments with respect to the Nusselt number. Our results with respect 
to the r.m.s. temperature fluctuation field, 9 , are in only fair agreement; the 
theoretical estimate of < 9 2 > 1/2 is too large near the boundary layer by ~40%. 
This overestimate of 9 may be connected with the low level of turbulence activity 
predicted theoretically. It is not clear to what extent this particular error is due 
to the use of cr = °°, and how much is due to an inadequate numerical treatment 
of the problem. 

Our conclusions regarding the theoretical existence of an asymetric mean 
temperature field (of the sort found experimentally by Somerscales and Gazda) 
are not so simple. The calculations indicate that an initial asymetry of the mean 
field persists for a time long compared to a circulation time of a fluid particle. 
Nevertheless, there is a clear tendency in the present calculation for the asymetry 
to disappear, or at least to be reduced to a level much below that found by the 
above experimenters. For those runs for which the initial data does not induce 
any initial asymetry in the mean temperature field no measurable asymetries 
developed. Such cases include those runs started near the conduction state with 
completely random initial data. These results indicate that the hexagonal type 
plan form is not stable (against finite perturbation). The predicted plan form 
(insofar as one can speak of such) appears to be one in which the motion is 
three dimensional, and consists of two interacting rolls, and suggests that found 
by Rossby for large a fluids if R > 10 4 . The shape of the mean temperature 



field is similar to that predicted by the quasilinear theory. The results at large 
R (R > 0.5 x 10 6 ) have the overshoot region just exterior to the boundary layer. 
However, the region is more diffuse in the present calculations that in the mean 
field calculation. 

H. THEORY AND INTEGRATION PROCEDURE 

(a) The equations of motion and their Fourier decomposition . 

The Boussinesq-Navier-Stokes equations at infinite Prandtl number may be 
reduced to the following nondimensional form: 


V • v = 0 


(1) 


-V 4 V = R 



(2) 



<9 = /3w - V • ( v (9 - v i9) 


(3) 



(4) 


Here, v in the fluid velocity field, T = 6 +T is the temperature field, T is the 
horizontal average of the temperature field, R is the Rayleigh number, k is a 
unit vector in the vertical direction, /3 is (-dT/dz), and w is v z . The infinite 
conduction - rigid boundary conditions imply that (Chandrasekhar (1960) ): 
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T(0, t) = 0 


T(Z = 1, t) = -1 


w(0, t) = g (0, t) = w(l, t) - ||(l,t)=0. 

These equations for T are appropriate for perfectly conducting boundaries. The 
nondimensionalization is such that the top plate is at Z = 1. We do not impose (ex- 
plicitly) any lateral boundaries, in view of the fact that we are interested in the 
solution of (1) - (4) under the condition of statistical (horizontal) homogenuity 
and isotropy. A convenient way of treating such a system is to use Fourier 
transform techniques, in the horizontal. We use it also in the vertical. 

Let 6 “ denote the component of Q (r) proportional to sin n^Z e rT7a ’ (l x + ie), 
and let /J n be the component of -(dT/dZ) proportional to cos rm Z. Then 6* and 
P n satisfy, 
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A derivation and discussion of these equations have been presented elsewhere 
(Herring (1964) ). Here we only note that the J n a m factor serves to eliminate the 
velocity field from the cr = °° system. The last term in the definition of J repre- 
sents the effect of the rigid boundary conditions. 

(b) Initial Value Problem for Homogenious and Isotropic System 
Consider those solutions of (1), (2), (3), and (4) which possess horizontal 
homogeniety and isotropy. By this we mean that the spacial average of the v - T 
fields (averaged over a finite but large area A) are independent of the horizontal 
position of A, and depend only on the vertical distance from the lower boundary. 
The linear dimensions of A should be large compared to correlation lengths such 
as 

l = J T(x) T(x + L) dL/J(T(x)) 2 dL. 
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To set up such solutions as an initial value problem we employ the Fourier am- 
plitude equations for 0 n a and T n . First the wave number space a is descritized 


to a i, and (5) and (6) are replaced by a finite set of N equations for 0“ 1 and T n . 
At t = 0, the 0^ 1 are assigned random values by selecting each complex number 
0“* from a Gaussian set of random numbers. The randomization is with respect 
to ai only and correlations between modes 0* 0“ a , (n ^ m) are permitted. In- 
deed, such correlations are necessary if the initial temperature field boundary 
layers. If we now take the limit of such a system as N — 00 (in such a way that 
the a F s "fill up" the entire finite wave number plane) then the initial tempera- 
ture field T (r, 0) becomes spacially homogenious and isotropic. We note that 
this method of posing the problem is one for which horizontal averaging is iden- 
tical to ensemble averaging over the initial set of randomly selected amplitudes. 

In practice we must always deal with a finite number of the ai_'s, and rFs 
(vertical wavenumber). With the above procedure for selecting initial data, the 
truncated system can only achieve approximate isotropy and homogeniety. More- 
over the subsequent development of the system will be spuriously affected by the 
finite dimensionally of the system (N). Suppose, for example, a particular wave 
number in a finite system is (accidentally) strongly excited at the initial time. 
This mode (or perhaps another to which this mode transfers its energy by non- 
linear interactions) will eventually have a dominant role in the system. If the 
number of degrees of freedom of the system were very large, then the chance 
over excitation of any mode would be correspondingly small; consequently it 
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would take a very long time for it or the mode to which it transfers its energy it 
feeds to become dominant. 

In the mean time the average properties of the flow come to quasi-steady 
equilibrium. It is the properties of the system during this quasi-steady state 
which interests us. In the limit as the number of modes becomes infinite, this 
quasi-steady state becomes a time stationary state. No mention has been made 
here as to whether the flow is turbulent. Presumably, if the experimental flow 
is turbulent the above method will give the turbulent results. If however the flow 
is cellular and static, the present calculation may differ in principle from ex- 
periments. 

Because of practical limitations in the size of the calculations, it is impor- 
tant to have some approximate measure of the degree of horizontal isotropy of 
the T-field. It is measured by a two dimensional tensor, 

'■‘ = f dn ik i 


Here, 


Vj = i 3/Bi + j B/BXj , and Xi , Xj 


are horizontal orthogonal coordinates. The volume of integration is over the 
entire field. An elipse may be associated with I in the usual way. If the flow is 
isotropic in the Xj, X 2 plane, the semimajor and semiminor axes are equal, for 
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flows which only depends on one dimension, the semiminor axis is zero. We 
define an isotropy parameter p as 

p = 1 - (a - b)/(a + b) . 

Here a is the semimajor axis and b is the semiminor axis. The p-parameter 
does not completely characterize the angular distribution of the flow, since p = 0 
for the nonisotropic hexagonal or square plan form also, p characterizes the 
second moment of the angular distribution only. A complete characterization of 
the angular distribution, needs a specification of the rest of the moments. 

(c) Method of Truncating the Wave Number Spectrum . 

The truncation procedure used here to reduce the infinite system of equations 
for 0* and fi n to a finite set guarantees two important properties of the flow: 

(1) entropy conservation and (2) the prescription must be appropriate to describe 
an approximate horizontally isotropic system. Entropy conservation is 
simplest to guarantee: one may verify from equations (5) and (6) that any trun- 
cation prescription which selects an arbitrary set of a i and the first N vertical 
wave number for and the first 2N vertical wave number for fi n satisfies a 
finite wave number version of the entropy conservation law: 


_d_ 

dt 



2N N 2N 2N 

2^iTi 2 +2Z v " |s ;i 2 +Z>" IT " |2 = 1+ ”I] nT » 

n n , a n= 1 n= 1 


( 7 ) 
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In applying the above truncation prescription, the convolution term in (5) is taken 
to be non zero only if the three wave numbers a, a , a" satisfy a = a' + a". 

How to choose the set ai so as to best achieve the goal of isotropy in the hori- 
zontal is not so clear. Two possible methods suggest themselves. 

The first method is prescribed by generating a set of a l Is by the equation 

Sj = ( i n + jm) C, (n, m) = i = (0, 1, 2, — ) « (g) 

and then discarding any a i for which 

|ctj_ | < a Q , or |aj | > a r 

Here a Q and a l are arbitrary cut off wave numbers and c is an arbitrary constant. 

For a calculation to be realistic a Q and must be outside the entropy containing 
region of the spectrum. It is clear that for a finite system the above procedure 
favor square plan forms. 

A second method consists in first generating a set of constant length vector 
by means of the formula, 

ap = (i cos — + i sin a. , p = 0, 1, • • • N* 

and then constructing all sums and differences a q of this set for which a Q < |aj < a r 
This gives exact isotropy under a finite rotation of angle (P/N*)7 t, P = 0, 1, 2, ... N* , 
and does not favor any plan form unless N* is divisible by 4 or 6, for which 
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square or triangular (hexagonal) plan forms are preferred. This second method 
gives a larger number of nonlinear terms than the first for a given total number 
of modes included in the system. 

(d) Numerical Integration Procedure 

The numerical integration method consists in a modified second order pre- 
dictor corrector technique. This relatively low order of the scheme is dictated 
by machine storage considerations. The modification referred to consists chiefly 
of a stabilization of the high wave number modes. Equations (5) and (6) may be 
symbolically written as 

Ji + » X •= F(x). 

A formal integration of this equation yields 

X(t + A) = X(t) e" vA + f dt' e~ vt ' F(X(t')>. 

n> 

From this equation we extract the following (constant A) second order predictor 
correction scheme: 


y (t + A) = X(t) e- vn + 


vA , (l-e- vA ) 


2v 


{ 


(l-a)F(X(t)-A))+ aF(X(t)) 


( 10 ) 


X(t + A) = X(t) e- vA + £! 6 


£4 


F(X(t)) + F(y (t + A)) 


( 11 ) 


13 




Here y ( t + A) is the predicted value of X ( t + A) ; a is an as yet arbitrary num- 
ber, whose value affects the stability (but not the order of accuracy) of the scheme. 
If a = 1 the above procedure (except for the e" v ^ -factor) is the standard Adams- 
Moulton method. 

Methods based on (10) and (11) are not entirely compatible with the entropy 
constraint (7). To see this, consider the simple case in which v = 0, and T = 0. 
Then (7) implies that 


|3 a i 2 

1 n 1 

n , a 

is a constant of motion, while (10) and (11) give 

2 \e a n (t + A) I 2 = 2 l0£(t)| 2 + 0(A 3 X). 

We chose a so that the long time average of the entropy constraint is most nearly 
satisfied. Clearly one way of doing this is to pick a value of a each time step 
which forces (7) to be satisfied. However, this is too complicated to be practical, 
and some numerical experimentation persuaded us that a = -1/2 is the best, 
simplest overall choice. Nordsieck has pointed out that the choice a = -1/2 (for 
second order integration schemes) best preserves the finite time step reversi- 
bility condition for reversible systems. He also points out that this scheme is 
the most unstable. However, in our system, the present scheme is stabilized by 
conductivity factors e“ v ^ serve to stabilize. 
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An integration technique which preserves the entropy constraints exactly is 
to transform dependent variables from the (complex) <9“ field to variables E 
and 0, where 6 = Ve" e 11 ^. Such a scheme was tried, but gave excessively large 
error for E near zero because of stability problems. 

IV. RESULTS AND DISCUSSION 

(a) Mean Properties of the Flow 

Results for the mean temperature field T(Z) and the r.m.s. temperature 
field 0(Z) = (<T 2 (r,t)> - <T> 2 ) 1/2 are presented in figures (1) - (2) for a range 
of Rayleigh numbers R = 4 x 10 3 , 10 4 , 10 s , and 7.39 x 10 s . These results are 
both horizontally and time averaged. The horizontal wave number spectrum is 
determined by the second method described in section IIIC, and the parameters 
are, N* = 6, a Q = 1 J/2, aq = J/2\ The number of vertical wave numbers 
ranges from 5 forR = 4000 to 20 for R = 10 6 . We postpone till section IVC a 
discussion and estimation of the wave number truncation error. The time 
averaging of the 0 and T fields was carried out only over the last portion of the 
time displacement of the system - after the excursion in N u had become less 
than 10%. Time integrations were stopped when the time average of N u become 
steady to three significant figures. The isotropy parameters for the R = 4 x 10 3 
and R = 10 s runs are shown in figure 3. 

Results for the averaged Nusselt number are shown in figure 4. The range 
of Rayleigh number covered does not warrant an extraction of a power law from 
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the data. However, if one assumes a N u ~ R 1/3 law (as has been analytically 
forecast by Roberson for two dimensional convection) then a continued fraction 
extrapolation of the data predicts 

N _ 0.945 R 1/3 for R co . 

U 

The points in figure 4 give the experimental results of the author cited in the 
introduction. 

(b) Comparison of Present Calculation with Quasilinear Theory and with 
Experiments 

It is of interest to compare the present results with those of the quasilinear 
theory, which simplifies the basic equations (1), (2), and (3) by deleting the 
last term in equation (3). A comparison of the two calculations is presented 
in figure (5) for the mean temperature field T(z), at R = 10 5 . For the 
quasilinear system, N y = 11.5 while for the complete system N u = 9.18. The 
main difference between the two calculations is seen to be that the overshoot 
region, just exterior to the boundary layer is much more diffuse according to the 
present calculation than according to the quasilinear calculation. The inclusion 
of the fluctuating-self interactions does not completely remove the overshoot 
region, as was one time thought (Herring (1964)). 

The numerical results may also be compared to experiments at high Prandtl 
number. Figures 6 and 7 give a comparison of the present results for T and 
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i p l/2 at R = 7.39 x 10 5 with the recent experimental results of Somerscales and 
Gazda. The most important differences between the present results and the ex- 
perimental results is the fact that the present results for T and 0 are very 
nearly symmetric about the mid-plane of the flow. Such a high degree of sym- 
metry is not a result of time averaging but is also obtained at any instant pro- 
vided the system is near the steady state. This disparity is curious. Within 
the framework of the Boussinerg approximation, an asymmetric T and 0 fields 
may be expected if the large scale part of the flow field is nearly statis and if 
the plan form of the motion is more or less hexagonal. By this we mean that 
the 9(x, y, z)-field has its large scale (vertical) structure dominated by modes 
(P , d a , and 6 a " (a = a* + a") such that 9 a , 9 a ' and 6 a " are in phase (for an 
exact hexagonal plan form, Q a - 6 a ' = 9 a "). The asymmetry corresponds to 
ascending (descending) columns of fluid of high velocity but low crossectional 
area and descending (ascending) columns of fluid of small velocity but large 
crossectional area. The direction of asymmetry is arbitrary in the Boussinerg 
approximation and is determined by the slight dependence of viscosity on the 
density. Our calculations admit such a solution-provided these solutions are 
stable to finite perturbations. Apparently these solutions are not stable in as 
much as the solution we find passes a near symmetric T field. Our solutions 
compound to having the large scale part of the 9 ( x, y, z) field discribable by 
modes 6 a , 9 a ', and 9 a " (a = a' + a") such that 9 a ‘ ~ 9 a " , and 9 a being (ap- 
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proximately) the harmonic of Q a ‘ and d a " . If only these three modes are in- 
cluded in the calculation, and if the solution is static, then it may be verified 
that the solution has exact mid-plane symetry. This solution corresponds to 
two roll structures superimposed at an angle of 30°. A possible explanation of 
the disparity is that Somerscales and Gazda actually sampled the temperature 
along a line in a convective cell, which remained fixed in space during the 
course of the measurement. 

The computed 0 field shown in figure 7 overestimates the r.m.s. temperature 
fluctuation field near the boundary layer by about 40%. This overestimation is 
probably due to the absence of strong small scale turbulence in the numerical 
experiment, and to its presence in the real experiment. What is not clear is 
how much of the difference is due to the fact that our numerical results are at 
infinite a , and the experiment at finite a , and how much is due to the inadequate 
numerical representation of the small scale structures. 

The computed values of the Nusselt number are in good agreement with 
experiment (N v (R = 7.39 x 10 s ) = 9.03 compared to our 9.18). However, this 
agreement is somewhat accidental, and depends on our arbitrary choice of a Q . 

We shall examine in section I VC the dependence of the T -field on a Q ; it turns 
out that the value of N v = - (dT/dz ) 0 depends more sensitively on a Q than 
(dT/dz ) at other values of Z. 
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(c) Discussion of Accuracy of the Calculation 

The two important sources of numerical errors is the present calculation 
are time descritization error and wave number truncation errors. To eliminate 
the time step errors, the time step was decreased until the evolution of 
certain selected quantities whose accurate evolution was desired were predicted 
with tolerable accuracy. These quantities were the Nusselt number, the isotropy 
parameter, p, and the entropy spectrum defined by 


S 


a 

n 



e a e~ a . 


Naturally, this means that certain other quantities, for example 0“ for a and n 
in large, are not predicted with great accuracy. The time steps were, At = 0.01, 
for R = 4,000; At= 0.025, for R = 10 4 ; At = 0.0025 for R = 10 5 , and At = 0.001, 
for R = 7.39 x 10 s . 

Wave number truncation errors exist both with respect to the vertical cut 
of wave number, n 0 , and with respect to the horizontal wave number cut offs 
a Q and a . The appearance of significant vertical wave number truncation error 
is manifest by spurious oscilations in the mean field of frequency 7m 0 . Such 
errors thus easily recognized are easily eliminated. 

Assessment of horizontal wave number truncation error is more difficult, 
since such errors do not manifest themselves in terms of easily recognizable 
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unphysical behavior for 0(7, t) or T(7, t). For example, the 0 and T fields for 
the quasilinear system behave quite reasonably despite the fact that these ampli- 
tude contain only one a . It appears plausable that these errors are eliminated 
for a given vertical wave number n if S n ° and S n 1 are much smaller than 
the peak value of S a : in other words the spectrum S“ should be adequately 
represented within a wave number band a Q , . Data on this point is presented 
in Table I for R = 7.39 x 10 s . From these data, it appears that there is little 
truncation error near a Q . Also, there appears to be little error in the large 
scale part (the odd modes) of the s “-field. However, the even modes appear to 
have an appreciable error, since near a 1 these intensities are maximum. These 
errors, however produce little change in the mean field or 0 (z)-field. 

Another assessment of the a-truncation error may be obtained by rescaling 
the a -spectrum by multiplying each value of a by a constant factor. Such a 
scale change amounts to selecting a new (and in principle exact) solution to the 
T , 9 system. If both of the a spectra cover the energetically significant portion 
of a-span one expects little difference in the horizontally averaged properties 
of the flow. However, our results cover only a small annulus in a -space, so 
that it is important to have some measure of the difference. Figure 8 and 9 
give a comparison of the T and 0-fields for a Q = 1/ / 2 and a Q = 1.0. The 
respective values of N u are 9.18 and 10.21. We observe that the two mean fields 
are indistinguishable, with only slight differences in the 0 fields. Apparently N u 
is the most sensitive to a rescaling of the a spectrum of all the mean field 
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quantities. This fact seems to vitiate any claim to predicting accurately theN 


versus R curve at large R. However, we note that of T and i p are insensitive 
to the above variations of a Q . 

V. CONCLUDING COMMENTS 

In the present paper we have developed numerical solutions to the full 
scaled three dimensional thermal convection problem, at infinite Prandtl number, 
and for rigid boundaries. The procedure used was designed to give an approxi- 
mate account of horizontally homogeneous and isotropic flow situations. To 
implement this we used the Fourier transform of the equations of motion. Our 
results, which included a maximum of 36 horizontal wave number vectors and 
20 vertical wave numbers appear to adequately describe the flow fill up to 
R = 10 5 ; beyond this R the results appear to show horizontal wave number trim- 
cation error. This error seems to affect the boundary slope of the mean 
temperature field more than other mean quantities, such as T and i p. Despite 
some numerical uncertainties, certain of the qualitative features of the flow fill 
are predicted with reasonable confidence. 

The mean temperature field T developes a negative gradient just exterior 
to the boundary layer, similar to that predicted by the quasilinear method. Thus 
the presence of the fluctuating-self -interactions does not entirely remove the 
negative gradient region, but does weaken it and makes it more diffusive. This 
type of mean temperature field is probably characteristic of large a fluids, and 
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appears to result from the fact that at large a the flow is characterized by large 
scale velocity field sweeping small scale temperature structures into the 
boundary layer. 

We do not find the large asymetrics in the mean temperature field reported 
by Somerscales and Gazda. Moreover, it does not appear likely that this result 
depends on wave number truncation error, since this behavior is primarily con- 
nected with the large scale part of the T-field. Possible reasons for this dis- 
crepancy are (1) the possibility that Somerscales and Gazda measured not the 
horizontal mean temperature field but rather the temperature field at a fixed 
vertical position in a slowly evolving convective cell, (2) in the experiment, 
the vertical boundaries play a significant role, (3) the finite (but small) 
value of Reynolds number in the experiment affects the results. At any rate, it 
appears that the hexagonal type asymetry in the mean field is definitely not 
stable with respect to finite perturbations. 

The dynamic of the flow as emerge from the numerical calculation consists 
of almost static large scale cells with more rapidly fluctuating small scale 
structures superimposed on them. The flow field is three dimensional (it does 

not degenerate into rolls) large scale part converts (roughly speaking) of rolls 
superimposed at an angle of 30°. The precise value of this angle depends on our 
method of selecting the a -spectrum, and no great significance shall be attacked 
to its value. 
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Table 1 


cl) 


n/a 0.7425 

0.3877E-04 

0.1707E-04 

0.6827E-05 

0.1734E-04 

0.6862E-04 

0.6419E-05 

0.1486E-04 

0.3482E-05 

0.2083E-05 

0.3422E-05 

0.1507E-05 

0.5182E-06 

0.5152E-06 

0.4823E-06 

0.1824E-06 

0.1501E-06 

0.5348E-07 

0.2157E-07 

0.4877E-07 

0.1041E-07 


1.0253 

0.7914E-04 

0.2105E-04 

0.1808E-03 

0.3224E-04 

0.1072E-03 

0.1641E-04 

0.1586E-04 

0.1072E-04 

0.8242E-05 

0.1008E-05 

0.2558E-05 

0.1501E-05 

0.8817E-06 

0.6162E-06 

0.5188E-06 

0.6226E-07 

0.1044E-07 

0.4123E-07 

0.1338E-07 

0.2876E-08 


1.2374 

0.3942E-04 

0.4508E-04 

0.5760E-04 

0.7151E-04 

0.8646E-04 

0.1684E-04 

0.3266E-04 

0.1386E-04 

0.4788E-05 

0.5775E-05 

0.2116E-05 

0.1066E-05 

0.7412E-06 

0.5470E-06 

0.2070E-06 

0.1203E-06 

0.6024E-08 

0.1550E-07 

0.2043E-08 

0.6730E-08 


1.3789 

0.1644E-04 

0.3409E-02 

0.5694E-04 

0.8746E-03 

0.1060E-03 

0.5760E-03 

0.2946E-04 

0.4128E-03 

0.4889E-05 

0.1768E-03 

0.3358E-05 

0.6728E-04 

0.5385E-06 

0.1336E-04 

0.1040E-06 

0.2557E-05 

0.1086E-06 

0.1700E-06 

0.2604E-07 

0.1404E-06 


Total 

0.1738E-03 

0.3492E-02 

0.3022E-03 

0.9957E-03 

0.3683E-03 

0.6156E-03 

0.9284E-04 

0.4409E-03 

0.2000E-04 

0.1870E-03 

0.9539E-05 

0.7036E-04 

0.2677E-05 

0.1501E-04 

0.1012E-05 

0.2890E-05 

0.1786E-06 

0.2483E-06 

0.9024E-07 

0.1604E-06 
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Figure 2. Mean squared temperature fluctuation as a function of z R - 4*10^, 10^, lO'’, 

and 7.39* 10 5 . 
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Figure 3- Isotropy parameter (see text), p as a function of time for R = 4*10 3 , and 10 5 . 
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Figure 5. Comparison of present mean temperature profile, T, with quasi-linear 
theory. Here R = 10^. 
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